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Abstract 

The evolutionary effects of glacial periods are poorly understood for Southern Hemisphere marine intertidal species, 
particularly obligatory sessile organisms. We examined this by assessing the phylogeographic patterns of the southern 
African volcano barnacle, Tetraclita serrata, a dominant species on rocky intertidal shores. Restricted gene flow in some 
geographical areas was hypothesized based on oceanic circulation patterns and known biogeographic regions. Barnacle 
population genetic structure was investigated using the mitochondrial cytochrome oxidase subunit 1 (COI) region for 410 
individuals sampled from 20 localities spanning the South African coast. The mtDNA data were augmented by generating 
nuclear internal transcribed spacer 1 (ITS1) sequences from a subset of samples. Phylogenetic and population genetic 
analyses of mitochondrial DNA data reveal two distinct clades with mostly sympatric distributions, whereas nuclear analyses 
reveal only a single lineage. Shallow, but significant structure (0.0041-0.0065, P<0.01) was detected for the mtDNA data set, 
with the south-west African region identified as harbouring the highest levels of genetic diversity. Gene flow analyses on 
the mtDNA data show that individuals sampled in south-western localities experience gene flow primarily in the direction of 
the Benguela Current, while south and eastern localities experience bi-directional gene flow, suggesting an influence of 
both the inshore currents and the offshore Agulhas Current in the larval distribution of T. serrata. The mtDNA haplotype 
network, Bayesian Skyline Plots, mismatch distributions and time since expansion indicate that T. serrata population 
numbers were not severely affected by the Last Glacial Maximum (LGM), unlike other southern African marine species. The 
processes resulting in the two morphologically cryptic mtDNA lineages may be the result of a recent historical allopatric 
event followed by secondary contact or could reflect selective pressures due to differing environmental conditions. 
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Introduction 

Past climatic changes in the form of glaciations, as well as 
associated changes in temperature and precipitation, have been 
suggested as drivers of demographic change in both terrestrial and 
marine species. Numerous marine species have been shown to be 
affected by sea level fluctuations that occurred during the last 100 
000 years and especially around the Last Glacial Maximum 
(LGM, 26 000 to 19 000 years ago) [1], [2]. Documented changes 
have included population bottlenecks and alterations of gene flow, 
as well as extinctions of local populations [3], [4], [5], [6], [7], [8]. 
In contrast, refugia provided a means of persistence for some 
populations, especially in species in the northern hemisphere [6], 
[9], [10], [11], [12], [13], [14]. The effects of sea level changes 
become less obvious in areas that had no ice cover during 
glaciation periods, such as southern Africa. Although the 
documented drop of — 1 30 m in sea level caused large areas of 
the Agulhas Bank to emerge with the formation of the southern 
coastal plain [2], [15], the effects on the topology of the southern 
African coastline are not well known. Further, no evidence for 
vicariance events, such as those shown along the south-western 
coast of Australia as a result of the land bridge across the Bassian 



Isthmus, [16], [17], [18], have been documented [19]. However, a 
number of southern African marine species show some degree of 
population genetic structuring along their geographic ranges, [20], 
[21], [22], [23], [24], yet barriers vary between species and there 
are few congruent patterns [19], [25], even in closely related 
species [26]. The variability in these patterns makes the 
explanation of the processes involved in generating the population 
genetic structuring particularly complex [26] . 

Such variance in patterns also applies to the study of 
demographic change; some marine species have been shown to 
maintain demographically stable populations throughout the 
LGM, without significant changes in population sizes, while 
others show population bottlenecks with post-LGM recoloniza- 
tion; both scenarios seem equally plausible [8] and can be found 
even in closely related species [26]. To better understand the 
evolutionary processes that drive changes in population sizes and 
distributions, it is important to carefully consider the timing of 
population expansions and contractions [27]. Further, it has been 
hypothesized that species with different life-histories show different 
patterns in response to climatic oscillations [8], [18]. For example, 
sessile organisms such as suspension feeders may be more likely to 
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persist in cooler climates due to the greater physiological costs and 
potential harm of mobile lifestyles under harsh conditions [8] . In 
addition, it has been suggested that the majority of species showing 
results consistent with LGM persistence are found in the lowest 
region of the intertidal or shallow subtidal and possess planktonic 
larvae [8], [18]. 

As with a number of studies, the link between life-history 
characteristics and population genetic structure is tenuous at best. 
Evidence from direct and indirect estimates of dispersal suggest 
that larvae consistently disperse on smaller spatial scales than 
expected based on their pelagic larval duration (PLD) [28], [29], 
[30]. Species with planktonic larval stages do not always result in 
panmictic populations [21], [31], [32], [33], [34], and species with 
poor dispersal abilities, such as direct developers or species that 
occupy small home ranges can show connectivity [16], [35]. 
Furthermore, a long PLD can contribute towards genetic 
panmixia [35], [36], [37], [38], [39], [40], [41], but known 
barriers to gene flow, which may be linked to oceanic circulation 
in a region and or to historical processes [3], [22], [33], [36], [42], 
[43], [44], can also play significant roles in shaping genetic 
patterns. 

In an attempt to contribute to the growing body of literature 
describing genetic patterns in southern African marine species and 
to determine the effect of the LGM on an intertidal organism, we 
chose the volcano barnacle, Tetraclita serrata, one of the 
dominant invertebrates of the mid intertidal zone with a wide 
distribution ranging from KwaZulu-Natal (east coast of South 
Africa), extending into Namibia (west coast of southern Africa; 
Fig. 1). Although this range spans at least four of the five 
biogeographic provinces described for the region, the species is 
more commonly found on the south and east coasts of South 
Africa (Fig. 1) [45]. A systematic study, incorporating both 
molecular and morphological data, suggests that T. serrata consists 
of two sympatrically distributed populations in South Africa [46]. 
Some conclusions were made on the demography of the species 
based on molecular diversity indices, but in the absence of 
molecular dating, the effect of the LGM could not be established. 
Furthermore, although hypotheses on gene flow were presented, 
no gene flow analyses were conducted. 

The present study provides a more in-depth phylogeographic 
consideration, based on a considerably larger data set, to elucidate 
the intraspecific variation of T. serrata across sampling localities 
along the South African coastiine. As barnacles are obligate sessile 
organisms with planktotrophic larvae, it is hypothesized that the 
longevity and dispersal capabilities of juvenile stages and the 
influence of the associated current systems may dictate the 
dispersal and genetic connectivity of the species. Given their 
position in the mid-intertidal, it is also postulated that T. serrata 
was not severely affected by the LGM (see [8], [18]). To test these 
hypotheses, the cytochrome oxidase subunit 1 (mtDNA COI) and 
the nuclear internal transcribed spacer 1 (nDNA ITS1) were used 
to document the phylogeographic structure of T. serrata; gene 
flow patterns and demographic changes of the species were 
investigated using the mtDNA data. 

Materials and Methods 

Study area 

The oceanographic patterns of southern Africa are complex and 
dynamic with seasonal effects on offshore, as well as coastal marine 
species [47]. The area is governed by two major boundary 
currents; the cold-water Benguela Current, which flows north- 
wards along the west coast of the continent, and the warm-water 
Agulhas Current, which moves southwards along the east coast of 



South Africa (Fig. 1). The diverse oceanographic and biological 
regimes along the South African coastline have been recognized to 
fall into five distinct biogeographic provinces, each characterized 
by a unique suite of potentially physiologically-adapted plant and 
animal assemblages [48]. 

Collection of material 

A total of 410 adult individuals of T. serrata were randomly 
collected across 20 intertidal sites (Fig. 1; Table 1). Sampling 
localities span most of the distribution and have been partitioned 
according to coastal region (west, south-west, south and east; 
Table 1), which broadly correspond to four of the five coastal 
bioregions as defined by Griffiths el al. [48]. All samples were 
collected in the form of whole individuals and stored in 100% 
ethanol until DNA extraction. Sampling was carried out under 
permits from the Department of Agriculture, Forestry and 
Fisheries (RES20 10/01, RES20 11/12). 

DNA extraction, PCR amplification and sequencing 

Genomic DNA was extracted from cirral tissue using the 
Wizard SV genomic DNA purification system (Promega). PCR 
amplification and sequencing of a fragment of the COI gene was 
initially carried out using universal primers [49]. The data 
generated by these primers were used to design species-specific 
degenerate primers which were subsequently used in all reactions 
(COI-TsF: 5'-CGG RGC TTG ATC CGC YAT RGT MGG-3' 
and COI-TsR 5'-GCT CCG GCT AAA ACT GGA AKH G-3'). 
Nuclear DNA data were generated from a subset of samples by 
amplifying the ITS1 region (including partial sequences of the 18S 
and 5.8S genes that flank the ITS1 region) using the primers of 
Bass et al. [50]. Amplification was performed in 50-ul reactions 
containing 4 Jul of 10-100 ng DNA, 2.5 (U.1 each primer (0. 1 mM), 
5 ul lOxPCR reaction buffer (Super-Therm), 5 ul dNTPs 
(0.2 mM, Thermo Scientific), 4 ul Mg 2+ (Super-Therm), 0.2 ul 
Taq polymerase (Super-Therm) and 13.4 ul distilled water. PCR 
reactions were carried out in an ABI GeneAmp 2700 automated 
thermocycler with the following conditions: an initial denaturation 
step at 94°C for 3 min, followed by 35 cycles of denaturing at 
94°C for 30 s, annealing (for 30 s) and an extension at 72°C for 
45 s, ending with a final extension at 72°C for 5 min. PCR 
annealing temperatures were set at 58°C and 56°C for the COI 
and ITS1 loci, respectively. The amplified products were gel 
separated and purified using the Wizard SV Gel PCR Clean-Up 
System (Promega). Purified products were cycle-sequenced using 
BigDye (Applied Biosystems [ABI]) and analyzed on a 3730 
automated sequencer (ABI). 

Genetic diversity 

Sequences were manually inspected and edited using BioEdit 
v7. 0.9.1 [51], and the mtDNA nucleotides were translated into 
protein codons to confirm functionality. All haplotypes have been 
deposited in GenBank with the following accession numbers 
KC935448-KC935857 (COI) and KC967663-KC967682 (ITS1). 

Haplotype designation of the nuclear DNA data was performed 
in PHASE, using a coalescent-based Bayesian method [52] as 
implemented in DnaSP v5.1 [53]. A total of 1 000 000 iterations 
were performed using default settings. Arlequin 3.5 [54] was used 
to calculate the number of polymorphic sites and the haplotype (h) 
and nucleotide (71) diversities for both marker types. 

Analyses of population structure 

Population genetic structure was investigated for the COI data 
using analyses of molecular variance (AMOVA) in Arlequin 3.5 
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Figure 1. Map showing the distribution of Tetraclita serrata from KwaZulu-Natal to Namibia. Sampled localities (see Table 1) for Tetraclita 
serrata are shown in addition to coastal regions [west (W), southwest (SW), south (S), east (E)] in which samples were collected. These coastal regions 
correspond broadly with the Namaqua, south western Cape, Agulhas and subtropical Natal Bioregions, respectively [47], The frequency of Clade A 
and Clade B within each sampling locality are indicated using pie charts. The major ocean currents are shown. The direction and intensity of gene 
flow between sampled localities, and the overall direction of gene flow are indicated by arrows; mutation-scaled migration rates between localities 
are indicated above these arrows. 
doi:1 0.1 371 /journal.pone.01 021 1 5.g001 



[54]. Differentiation across sampling localities was examined 
through pairwise (Dst comparisons. Bayesian geographical struc- 
turing was investigated using BAPS v5.4 [55], where the 
geographic coordinates of each sampling locality were obtained 
from Google Earth. To investigate geographic distance as a 
contributor towards genetic structure, isolation by distance [56] 
was tested using a Mantel test [57] as implemented in Arlequin 3.5 
with 1000 permutations. The geographical distances between 
sampling localities were measured as the shortest land-sea interface 
in Google Earth. 

Statistical parsimony haplotype networks were constructed 
using TCS 1.21 [58] and haplotypes were connected using a 
95% confidence interval. The higher level clustering of clades that 
were not connected in the haplotype networks was depicted using 
Bayesian analyses, coupled with the Markov Chain Monte Carlo 
(BMCMC) inference, conducted in MrBayes v3.04b [59]. The 
most likely (best probabilistic) model of sequence was determined 
using the Bayesian Information Criteria (BIC) and Akaike 
Information Criteria (AIC) for the COI and ITS1 gene, 
respectively, as employed in the program jModelTest vO.1.1 
[60]. Bayesian analyses of the COI and ITS1 sequence data were 
run for 1 0 million and five million generations, respectively, with a 



burn-in of 25%. For each gene, two different runs were conducted 
to test for convergence, each with four chains and sampling 
frequency of 1000. Thereafter, the post-burn in samples of the two 
runs were combined and a 50% majority-rule consensus tree was 
constructed in MrBayes v3.04b. Clade support was assessed using 
posterior probabilities for the Bayesian tree and only values &0.95 
were regarded as robust support. Tetraclita squamosa and 
Tetraclita pacifica were used as outgroup taxa (GenBank accession 
numbers DQ647770 and DQJ563746, and DQ363680 and 
DQ363740, respectively). 

Demographic history 

Arlequin 3.5 was used to test the influence of historical sea level 
changes on the demography of T. serrata. Fu's F s test [61] was 
performed to test for deviations from neutrality. To confirm that 
departures from neutrality were caused by demography and not by 
natural selection, the McDonald-Kreitman (MK) test was 
performed in DnaSP v5.1, using T. squamosa and T. pacifica 
(GenBank accession numbers DQ363697-DQ363706 and 
DQ363686-DQ363695) as outgroups, respectively. This was 
followed by a mismatch distribution to test for population 
expansion and to estimate parameters needed to date possible 
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population expansions [54], [62]. The formula T = x/2 \xk [62] 
was used for dating approximations. As the exact mutation rates 
and generation times for T. serrata are unknown, we used the 
generation times of various barnacle species, ranging from one to 
two years [8], [32], [63], [64], and two published COI lineage 
mutation rates for barnacles (1.55xl0~ 8 [65] and 2.76xl0~ 8 
substitutions per site per generation [18]). 

A Bayesian Skyline Plot (BSP) [66] was constructed using 
BEAST vl.8 [67] to estimate changes in population size over time 
and to date a possible population expansion. The unique 
haplotypes of each COI data set (Glades A and B) were run for 
300 million and 200 million generations, respectively, under a 
GTR+I+G nucleotide substitution model with individual param- 
eters estimated from the data and a constant skyline model with 
five groups. For each BSP, prior distributions for the root height of 
the population were notified by initial estimates from the mutation 
rates used in the mismatch distribution analyses. The chain was 
sampled every 100 000 generations and the first 10% was 
discarded as burn-in. Trace plots were inspected to assess 
convergence, mixing and stationarity of the MCMC process in 
Tracer vl.6 [68]. The effective sample sizes were checked and 
confirmed to be &200 in order to avoid autocorrelation of 
parameter sampling. 

Connectivity and direction of relative gene flow between 
sampling localities 

The program Migrate 3.1.6 [69] was used to calculate the 
relative migration rates between sampling localities inferred from 
the mtDNA data. Given the linear South African coastiine, a 
stepping-stone model was used [70]. The settings were as follows: 
10 short-chains, each with a total of 2,500 generations and a 
sampling increment of 20 generations, and 3 long-chains, each 
with a total of 50,000 generations and a sampling increment of 50 
generations. The first 10,000 genealogies were discarded (burn-in). 
An adaptive heating scheme with four chains (starting values of 
1.00, 1.50, 3.00 and 6.00) and a swapping interval of one was used 
to ensure that efficient mixing occurred. For the other settings, 
default values were implemented. The analyses were run twice to 
ensure congruence. 

Results 

Genetic diversity and population structure 

Mitochondrial DNA analyses. The COI data analyses were 
based on a final alignment of 499 bp for 410 individuals. The 
maximum parsimony network, Bayesian tree and spatial clustering 
of all individuals (BAPS) indicate that T. serrata comprises two 
genetically distinct groups with strong nodal support. The two 
genetic assemblages occur sympatrically (Fig. 2a-c). These 
haplogroups represent two reciprocally monophyletic lineages, 
both showing Bayesian Posterior Probabilities of 1.00 (BPP). These 
are labeled Clades A and B (Fig. 2) and are separated by a mtDNA 
sequence divergence of ~8%. In proportion to Clade A, the 
frequency of individuals belonging to Clade B is markedly smaller 
(a sampling ratio of 12:1, Fig. 1). 

Within Clade A, haplotype diversity was 0.979 (±0.003), with 
187 haplotypes detected, and within Clade B it was 0.976 
(±0.019), with 26 haplotypes detected. The nucleotide diversity of 
the complete data set was 0.019 (±0.01); generally mtDNA genetic 
diversity was highest on the south-west and south coasts, with a 
progressive decrease towards the east coast and the lowest values 
on the west coast (Fig. 3a). With the removal of Clade B, the 
nucleotide diversity values are lower, with the south-west coastal 
region generally being the highest (Fig. 3). Although the haplotype 



diversities determined for each locality for the entire data set 
(Clade A and Clade B individuals) were generally high (>0.94, 
Fig. 3b), on a broader scale, sampled localities on the south-west 
coast, as well as the closest localities sampled on either side of the 
south-west coast region, namely Jacob's Bay, Kommetjie and 
Cape Agulhas, each had more than 50% unique haplotypes 
(Table 1). In contrast, those haplotypes shared among regions 
occur more frequently along the east coast than in the other 
regions (i.e. fewer private haplotypes, less than 45%). 

Within Clade A, the haplotype network (Fig. 2b) shows that a 
considerable amount of haplotype sharing is found along the east 
coast. Also, west and south-west coast sampled localities group 
together and are more divergent, suggesting spatial genetic 
structure. In Clade B, the majority of haplotypes are representative 
of individuals sampled along the west, south-west and south coasts, 
with only two Clade B individuals sampled as far as Haga Haga. 
Many of the haplotypes within this haplogroup are separated from 
the common haplotype by an intermediate "unsampled or extinct" 
haplotype (Fig. 2b). 

AMOVA suggests shallow, but significant, genetic differentia- 
tion among the 20 sampled localities (<J>st = 0.035, P<0.01); most 
of the variation was present within sampled localities (96.44% for 
the entire data set). When sampling localities were partitioned into 
their respective coastal regions (Table 1), AMOVA revealed 
significant differentiation among groups (C>ct = 0.041, P<0.01). 
When the individuals belonging to Clade B (Fig. 2b) were removed 
from subsequent analyses, the 0 ST value increased (0.065, P< 
0.01). For the entire data set, pairwise fl>sT analysis showed that 
Kommetjie was significandy differentiated from all localities east of 
Cape Agulhas (south coast), with 0 ST values ranging from 0.09 to 
0.23 (south coast, Table 2). Kommetjie was also found to be 
significandy differentiated from most of the localities sampled on 
the west coast (Table 2). Betty's Bay was found to be significandy 
differentiated from all localities east of Herold's Bay (south coast), 
with <1>st values ranging from 0. 1 to 0. 19 (Table 2). Within Clade 
A, Kommetjie and Betty's Bay were significandy differentiated 
from all sampling localities, except from each other and from Port 
Nolloth, with 0 ST values ranging from 0.07 to 0.34 (Table 2). Due 
to the relatively small sample size represented by individuals of 
Clade B, pairwise A>st analysis was not conducted for this clade. 

Nuclear DNA analyses. The aligned sequence data for ITS1 
were 426 bp in length, with 31 unique alleles and four indels 
detected. The ITS1 data displayed a high haplotypic diversity of 
0.960 (±0.019) and a low nucleotide diversity of 0.010 (±0.005), 
which is comparable with that of the mtDNA data. No significant 
structure was obtained and Clades A and B formed by the 
mitochondrial DNA analyses collapsed into a single haplogroup 
(Fig 2d, e). 

Isolation by distance 

A significant correlation between geographic distances and 
genetic differentiation was obtained for the entire T. serrata data 
set (Mantel test, Z = 8.925, P<0.01), revealing a larger Z-value of 
14.176 (P>0.05), which was not significant, with the removal of 
Clade B. The Mantel test was found to be non-significant when 
isolation by distance (IBD) was tested within each coastal region 
[west coast, Z = 0.119 (P>0.05); south coast, Z = -0.035 (P> 
0.05); east coast, Z = -0.096 (P>0.05)]. South-west coast localities 
were not included in the test because only two localities had been 
sampled in that region. 

Migration estimates between sampling localities 

Within the sampling range of T. serrata, the coalescent 
stepping-stone model revealed migration rates that were largely 
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Figure 2. MtCOl data as visualized in (a) a Bayesian tree (constructed in MrBayes MrBayes v3.04b), (b) haplotype network 
(constructed in TCS 1.21), and (c) spatial clustering of individuals (constructed in BAPS v5.4). NulTSI data is visualized in (d) a Bayesian 
tree and (e) a haplotype network. In b & e, the size and shading of the circles corresponds to the frequency of that haplotype and the coast from 
which the individual was sampled (west =W, south-west = SW, south =S, east =E). Each line represents one mutational step and perpendicular 
lines represent additional mutational steps. Each sampled haplotype is colour-coded according to the coast that the individual was sampled in. In c), 
each vertical bar represents one individual grouped by region and location in the same order as Table 1 . The shading of the vertical bars indicates to 
which group they are assigned by the BAPS analysis. 
doi:1 0.1 371 /journal.pone.01 021 1 5.g002 



bi-directional in the south and east coastal regions, whereas a more 
asymmetrical gene flow pattern was evident in a westerly direction 
between Mossel Bay and Hondeklip Bay (with the exception 
between Betty's Bay and Cape Agulhas) (Fig. 1). Interestingly, no 
southward gene flow was detected between Jacob's Bay and Betty's 
Bay, around Cape Point (Fig. 1). 



Demographic history 

For the mtDNA data, Fu's F s revealed negative values for the 
entire data set (although it was not significant for all sampling 
localities; Table 1), Clade A (F s = -17.858, P<0.01) and Clade B 
(F,= -25.065, P<0.01), respectively. The MK test was not 
significant for T. serrata when T. squamosa (G-value =0.841, 
P>0.05) was used as an outgroup. The more closely-related 
species T. pacifica could not be used for this test because of the 
lack of nonsynonymous fixed differences with T. serrata in the 
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Figure 3. Molecular diversity indices across the distribution range of T. serrata include (a) nucleotide diversity and (b) haplotype 
diversity values for the entire data set (Clade A+ Clade B, represented by diamonds) and for Clade A individuals (represented by 
squares). 

doi:1 0.1 371 /journal.pone.01 021 1 5.g003 



sequenced gene region. A pattern of population expansion was 
confirmed for both clades by the observed unimodal mismatch 
distributions (not shown), which were not significantly different 
from the Poisson distributions predicted by a sudden expansion 
model. This hypothesis of population expansion is supported by 
both the Sum of Squared Deviations (P>0.05) and Harpending's 
raggedness indices (P>0.05, Table 1), which tests for a fit to a 
unimodal mismatch distribution. 

The time since population expansion was calculated for mtDNA 
Clades A and B (Table 3). Analyses of Clade A suggest the time 
since expansion between ~40 000 and 145 000 years ago. 
Similarly, for Clade B, the time since expansion was calculated to 
be between ~46 000 and 164 000 years ago. Bayesian Skyline Plot 
analyses reveal dates similar to the mismatch distribution analyses. 
Analyses of Clade A show a population expansion took place 
between ~40 000 and 100 000 years ago, while analyses of Clade 
B show a population expansion took place between ~50 000 and 
100 000 years ago (Fig. 4). 

Discussion 

Within-population diversity 

Genetic diversity indices varied across the sampling range, with 
striking differences between sampling areas; generally, diversities 
were highest on the west and south-west coasts, and they decreased 
towards the eastern localities. A distinct peak in diversity was 
found for the south-west coast (Fig. 3). There were also notable 
differences in the distribution of private/ unique haplotypes, with a 
higher number of private haplotypes found in the south-west 
coastal region (Table 1). 

Regions of higher genetic diversity can be explained by several 
scenarios. Genetic variation can accumulate over evolutionary 
time if a region contains an area of high local retention or larval 
accumulation, especially where larvae accumulate in the near- 
shore environment [29], [71], [72]. High genetic diversity is also 
expected in a region where a large population has been 
maintained over evolutionary time, possibly due to the stability 
of the habitat [73]. If contemporary population size is an indicator 
of evolutionary history, then this will have contributed to the high 
genetic diversity as T. serrata usually maintains large population 
sizes on most rocky shores [45], [74], [75]. Irrespective of 
population size, regions of higher genetic diversity can also signal 
refugial areas during glacial low sea-level stands [76]. It can also be 
a signal of the absence of strong directional or stabilizing selection, 
as both mechanisms have been shown to reduce genetic diversity 
[77]. Populations at or near refugia should retain genetic diversity 
and may also be indicated by the presence of ancestral haplotypes 



[78]. Interestingly, Betty's Bay on south-west coast, possesses the 
highest combination of haplotype and nucleotide diversity 
(Tables 1, 2). A study on the commercially exploited rock lobster, 
Jasus lalandii, also identified the region between Betty's Bay and 
Kommetjie as one of high diversity [21], but the processes driving 
this area to be more genetically diverse than others are not 
understood. 

Genetic structure of Tetraclita serrata 

Population genetic structure of T. serrata is characterized by a 
significant AMOVA and a strong signal of isolation by distance for 
the mtDNA data. An important result was that T. serrata is 
comprised of two genetically distinct groups that are geographi- 
cally mixed (Fig. 2b, c), consistent with Tsang et al. [46]. Although 
Clade B is highly divergent and shares no haplotypes with Clade A 
at the mitochondrial level, the shallow population structure 
recovered reflects the geographic overlap of the haplotypes of 
the two clades, with the majority of the haplotypes being shared 
among the sampled areas (Fig. 2b). The pattern observed reflects a 
significant historical differentiation of the two clades and may 
represent an incipient speciation event. Lowered sea levels, as a 
consequence of historic glaciations, may have been responsible for 
the formation of allopatric lineages that have since come into 
secondary contact. 

In the absence of a clear allopatric, geographic explanation that 
could have given rise to the two clades in T. serrata, it is important 
to consider that physiological differences could also result in 
genetic diversification. Field observations and physiological 
experiments on the mussel Perna perna indicated that the western 
lineage of the species is intolerant of high water temperatures, thus 
explaining its exclusion from the east coast [79] . Similarly, in the 
mudprawn Upogebia africana, larvae of the eastern lineage are 
unable to tolerate cooler waters and, thus, are unable to establish 
themselves in cool-water provinces [80] . The mtDNA data of T. 
serrata suggest that Clade B extends no further east than the 
border of the sub-tropical and temperate regions at Haga Haga, 
whereas Clade A can be found in the sub-tropical region (Fig. 1, 
Fig. 2b). This implies that Clade B may be constrained to more 
temperate environments, as a result of physiological intolerance. 
The annual mean sea surface temperatures experienced along the 
South African coastiine ranges between 12°C on the west coast to 
26°C on the east coast [81]. Tetraclita serrata displays regular and 
annual recruitment and Dye [82] revealed a high variation in 
spatial and inter-annual recruitment intensity along the east coast 
of South Africa. It was also found that recruitment density was 
related to mean abundance [82]. This suggests that larval dispersal 
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Figure 4. Bayesian Skyline Plot of effective population size through time based on the mtDNA sequence data for Clade A (al and 
a2) and Clade B (bl and b2), using the mutation rates 3.1/My (equivalent to 1.55x10 8 substitutions per site per generation; al 
and bl) and 5.52/My (a2 and b2), respectively. The shaded area indicates the period of population expansion. 
doi:1 0.1 371 /journal.pone.01 021 1 5.g004 



and settlement may differ among regions. Along these lines, it is 
particularly notable that two distinct, species-level lineages of the 
Red Sea barnacle Tetraclita rufotincta occupy different vertical 
zones of the intertidal [83] . The genetic divergence in T. serrata 
may be explained by selective pressures due to ecological 
conditions, such as desiccation tolerance and predation. Given 
the geographic overlap of the two clades, it is possible that the 
genetic pattern found in T. serrata is the result of an ecological 



speciation process. Individuals belonging to Clades A and B may 
occupy different zones of the mid-intertidal (micro-habitat 
isolation), although this remains to be tested further. 

Despite the clade structure present in the mtDNA data and a 
sequence divergence of — 7.92%, the lack of geographic structure 
for Glades A and B and the overall lack of structure for the nuclear 
data suggest incipient speciation at best (Fig. 1, 2). This result is 
supported by Tsang et al. [46], who similarly found no genetic 



Table 3. The time since expansion for the entire mtDNA data set and mtDNA Clades A and B using expansion times (x) calculated 
in Arlequin 3.1 and published COI lineage mutation rates. 







Mutation rate 


Time since expansion 








k=1 


k = 2 


Entire data set 


1 = 3.457 


1.55x10 8 


111,740 


55,870 




2.76 x10~ 8 


62,753 


31,376 


Clade A 


t = 4.498 


1.55x10 8 


145,388 


72,694 




2.76 x10~ 8 


81,650 


40,825 


Clade B 


T = 5.078 


1.55x10 8 


164,135 


82,067 




2.76 x10~ 8 


92,178 


46,089 



Time since expansion and generation time (k) is represented in years. 
doi:1 0.1 371/joumal.pone.OI 021 1 5.t003 
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structure for nuclear DNA data (Histone 3) and no morphological 
differentiation. As the species is hermaphroditic, it can be ruled 
out that the lack of nuclear genetic structure is that of sex-biased 
dispersal, and is more than likely due to differences in the fixation 
rates of the two marker systems. 

Population genetic structure: associations with life- 
history characteristics 

Several phylogeographic studies for southern African marine 
species suggest the presence of at least four marine phylogeo- 
graphic breaks (see [19], [25]), of which we sampled across three 
(Cape Point, Cape Agulhas and the area between Port Elizabeth 
and the northern Transkei; see Fig. 1 of [19]). Although globally 
many studies have shown that life history characteristics such as 
pelagic larval duration (PLD) are not correlated with population 
genetic structure and gene flow [84], [85], planktotrophic larval 
dispersal of T. serrata may give some insight into the lack of 
structure recovered in this study. Even though the shared 
haplotypes found for the entire data set indicate that this brooding 
barnacle is able to disperse over long distances, suggestive of a long 
PLD, T. serrata does not form a panmictic population as several 
sampled localities are differentiated, even if <1>st values are low 
(Table 2). The presence of large numbers of haplotypes unique to 
certain regions (Table 1), particularly in the south-west coastal 
area (KO, WP, BB), provides evidence for population genetic 
structuring in the west and south-west coastal regions. The 
significant differentiation of Kommetjie (Cape Peninsula) and 
Betty's Bay (False Bay region) from other localities (Fig. 2b, 
Table 2) suggests that, at the finer geographic scale, dispersal may 
be restricted by water currents or by the local retention of larvae in 
the water column, as has been shown for other species [34], [86], 
[87], [88]. This is further supported by the signal of IBD for T. 
serrata. 

Gene flow in relation to ocean currents 

The strength and directionality of gene flow in T. serrata show 
interesting correlations with oceanographic features along the 
South African coast. All west coast (except Port Nolloth) sampled 
localities experience gene flow primarily in the direction of the 
Benguela Current, while half of the localities sampled on the 
south-west coast (see below) experience gene flow in the direction 
of the Agulhas Current (Fig. 1). Ocean currents and prevailing 
wind direction have long been recognized as important factors 
affecting the distribution, abundance and genetic variation of 
marine benthic invertebrates with planktonic larvae [89], [90]. A 
pattern of unidirectional gene flow, which is much lower than 
relative gene flow estimates around Cape Point, may play a role in 
the genetic structuring of sampled localities on the south-west 
coast, evident in the haplotype network with haplotypes from the 
west coast and south-west coast sampling localities clustering 
together (Fig. 1). Such reduced gene flow was also shown for rocky 
shore fishes [22], [26] and an endemic sea urchin, Parechinus 
angulosus [24] strongly suggesting that gene flow is reduced for 
inshore coastal species in this region. Cape Point is also 
characterized as a region where a major biogeographic break 
occurs [19]. The barrier is thought to be the result of the distinct 
temperature regimes on either side of the peninsula, which may 
influence larval dispersal and/or survival in the different 
environments. By contrast, the signal is more mixed on the south 
and east coast, with bi-directional gene flow, suggesting that both 
the inshore, wind-driven countercurrents and the offshore Agulhas 
Current play a role in the distribution of larval T. serrata [19]. 



Demographic history 

The combination of high haplotype diversity and low nucleotide 
diversity is a typical signature of population expansion following a 
bottleneck [91], [92]. Fu's F s simulations (Table 1) and the MK 
test, which was not significant, indicate that the observed 
deviations from equilibrium are not caused by natural selection 
and that T. serrata has undergone recent demographic change. 
The mtDNA haplotype network, showing several common 
haplotypes (Fig. 2), mismatch distributions (Table 1) and the 
timing of the most recent expansion event for T. serrata are 
indicative of a growing population. The population expansion of 
T. serrata predates the LGM, taking place 164 000 to 3 1 000 years 
ago, depending on the mutation rate and generation time used 
(Table 3). The BSP analyses agree with this (40 000 to 100 000 
years ago, Fig. 4). Although many southern African marine taxa 
show patterns of post-LGM population growth [8], [20], [21], 
[93], [94], [95], [96], also confirmed by a number of barnacle 
species found elsewhere [5], [8], several studies among rocky shore 
marine species also show pre-LGM population expansions (Pacific 
Ocean: [8], [18], [97], [98], [99], [100]; Atlantic Ocean: [4]). This 
result suggests that the dramatic and frequent climatic changes of 
the Pleistocene, pre-dating the LGM, may have had a stronger 
impact in shaping the demographic and biogeographic history of 
T. serrata [8], [101]. Even if a more conservative approach was 
used, such as the mitochondrial mutation rate based on the 
divergence rate of 2% per million years [102], the demographic 
changes of T. serrata would still have predated the LGM. The two 
mutation rates that were used from barnacle genera were 
substantially greater than the 2% rate and still the population 
expansion of T. serrata predated the LGM. Therefore, despite the 
limitations of calculating the timing of expansion [27], the 
variation in mutation rates alone cannot account for the relatively 
deep demographic history for T. serrata. 

Broadly, our results support Marko et al. [8] in that suspension 
feeders possessing planktonic larvae that occur in the mid- 
intertidal region show patterns of LGM persistence. The majority 
of studies to date on South African marine species showing 
patterns of post-LGM recolonization also possesses planktonic 
larval stages and is found subtidally, and is neither suspension 
feeders nor sessile species [7] . Nevertheless, these are certain traits 
worth considering for future research given that there are a 
number of barnacle species that possess planktonic larvae and live 
in rocky shore habitats showing clear population genetic signatures 
of regional persistence [4], [8], [18], [98]. 
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